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Abstract 

We consider the non-equilibrium dynamics of a chain of classical rotators coupled 
at its edges to an external reservoir at zero temperature. We find that the energy 
is released in a strongly discontinuous fashion, with sudden jumps alternated with 
long stretches during which dissipation is extremely weak. The jumps mark the 
disappearance of strongly localized structures, akin to the rotobreather solutions of 
the Hamiltonian model, which act as insulating boundaries of a hot central core. As 
a result of this complex kinetics, the ensemble-averaged energy follows a stretched 
exponential law until a residual pseudo-stationary state is attained, where the hot 
core has reduced to a single localized object. 

We give a statistical description of the relaxation pathway and connect it to the 
properties of return periods of rare events in correlated time series. This approach 
sheds some light into the microscopic mechanism underlying the slow dynamics of 
the system. 

Finally, we show that the stretched exponential law remains unaltered in the 
presence of isotopic disorder. 
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1 Introduction 



Many physical systems display unusual relaxation properties, such as slow 
kinetics in glasses [1] or dynamics of biomolecules [2], which are usually asso- 
ciated with the presence of a complex energy landscape. Long-living transient 
states also arise in non-linear systems, where the relaxation properties may 
be affected by the emergence of localized vibrations like solitons or kinks (see 
e.g. [3]). Of special interest is the role of a recently discovered class of spatially 
localized, time-periodic excitations, termed discrete breathers [4]. At variance 
with their counterparts in the continuum, they exist under very general con- 
ditions and their stability properties have been thoroughly investigated [4]. 
In particular, these objects may easily self-excite under very different phys- 
ical conditions, as exemplified by several numerical and experimental stud- 
ies [4,5,6,7,8,9,10]. 

A typical situation where breathers form spontaneously can be obtained by 
cooling the lattice boundaries, thus driving energy out of the system. In such 
conditions one is faced with a genuine nonlinear effect: the system prevents 
the complete dissipation of the energy by storing part of it in a residual long- 
living breather or multi-breather state. Such a phenomenon was first observed 
in a model of coupled harmonic oscillators pinned by an on-site nonlinear 
force [11,12] and later recovered in models of nonlinearly coupled oscillators 
without on-site interaction [13,14,15]. Nonetheless, the process of energy dis- 
sipation was found to exhibit quite different behaviors in these two classes of 
models. The main mechanism yielding a power-law decay of the energy to the 
residual state has been described in previous papers by the authors [13,14]. 
The stretched exponential behaviour, primarily observed in models with on- 
site interaction, was then conjectured to represent a typical feature for this 
class of models. Nonetheless, a satisfactory understanding of this peculiar de- 
cay process was still lacking. 

The present paper deals with the nonequilibrium transient dynamics of a 
model of coupled classical rotators on a one-dimensional lattice. Two dif- 
ferent cases will be considered, the one with pure nearest-neighbours coupling 
and the one with an additional on-site force which is a generalized version of 
the well-known discrete sine-Gordon chain. The quantum version of the lat- 
ter model has been introduced by Fillaux and collaborators [16] to describe 
inelastic-neutron-scattering experiments on methyl-pyridine crystals. As this 
system displays evidence of extremely slow kinetics [17], it is particully in- 
teresting to investigate how cooling and nonlinear effects may affect thermal- 
ization. Although limited to the classical case, our study may provide some 
insight on the observed phenomenology. 

The second motivation comes from the problem of stationary heat transport 
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in this class of models. It has been recently recognized that low-dimensional 
lattices show anomalous properties, namely that transport coefficients (e.g. 
the thermal conductivity) diverge [18]. The analysis of several models clarified 
that this striking feature occurs generically whenever momentum is conserved. 
For lattice models, this amounts to requiring that at least one acoustic phonon 
branch is present in the harmonic limit. Remarkably, the only known excep- 
tion is actually the pure nearest-neighbor model, which displays normal trans- 
port [19]. Since the transient dynamics is connected to the system's response, it 
is interesting to investigate the energy relaxation and to compare the peculiar- 
ities of each class of models. As we shall see, there are indeed some qualitative 
differences whose genuine nonlinear origin (excitation of long-lived localized 
structures) points at significantly differerent transport mechanism than those 
based on the customary phonon theory. 

The plan of the paper is as follows. The details of the system of coupled 
classical rotators cooled at its boundaries are presented in Section 2: we con- 
sider both the cases with and without on-site interaction. The simulations 
reported in Section 3 indicate that the energy relaxation dynamics displays 
in both cases a stretched-exponential behavior. The statistical interpretation 
presented in Section 4 relates such a stretched-exponential behavior with rare- 
event statistics. To our knowledge, this correspondence has never been pointed 
out before in the literature. In section 5 we show that the main features and 
the interpretation do not change when isotopic disorder is introduced. 

Two relevant results of our analysis merit to be stressed from the very begin- 
ning. The first one concerns the robustness of the phenomenon with respect 
to different variants of the model. The second one amounts to put such a phe- 
nomenon in a different perspective: the stretched-exponential energy decay 
process seems to be common to one-dimensional systems whenever the typ- 
ical solutions are static breathers, rather than being a peculiarity of models 
with on-site interaction. 



2 The model 

We study the relaxation toward equilibrium of a chain of classical rotators 
coupled on a lattice with damping acting at its edges. The equations of motion 
are 

lifa = -G sin fa + K [sin(0 m - fa) + sin(0;_i - fa)} - ~ffa [6 iA + 5 i:N ] (1) 

where fa in the angle variable of the i-th rotator (with % = 1, 2, . . . N). The 
Hamiltonian version of (1) is sometimes referred to as the sine-lattice equa- 
tions [20,21]. We impose free-end boundary conditions (fa = fa, (fijy = 4>n+i) 
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since we know from previous studies that this choice favors spontaneous local- 
ization [13]. In the following we will also consider the special case G = 0, where 
only interactions between nearest neighbors are present and refer to it as the 
NN model. Furthermore, we will set h = I for all i and K = 1.4911, a value 
that in our units corresponds to experimental data on 4-methyl-pyridine [16]. 

The general layout of a simulation has been detailed previously [13,14]. First, 
an equilibrium microstate is generated by letting the Hamiltonian (micro- 
canonical, i.e. 7 = 0) system evolve for a sufficiently long transient (typically 
500 time units). We used the 5 th -order symplectic Runge-Kutta-Nystrom al- 
gorithm [22]. The initial condition for the transient is assigned by setting all 
angles <pi to zero and by drawing velocities at random from a Gaussian distri- 
bution with standard deviation equal to the square root of twice the energy 
density E(0)/N. The resulting set of fa and fa is then used as initial condi- 
tion to integrate Eq. (1) with 7 > with a standard 4 th -order Runge-Kutta 
algorithm. The subsequent dissipative dynamics is described in the following 
section. 



3 The relaxation dynamics 

A typical simulation is illustrated in Fig. 1. In the left panels we show the 
space-time contour plots of the symmetrized site energies e,, defined as 

1 • 1 

e i = ^ + 2 K ^ ~ COS( ^ i + 1 ~~ ~~ C0S( ^ ~~ + G ^ ~ C ° S • ( 2 ) 

In the right panels we report the time series of the total normalized energy 
E(t)/E(0) for the same simulation runs. 

The energy density contour plot shows that the process is associated with the 
spontaneous creation of long-lived localized objects, characterized by a fast 
librating motion of a single rotator almost decoupled from its neighbors. These 
localized excitations (rotobreathers) have already been shown to exist in the 
Hamiltonian lattice [20]. The appearance of such strongly localized objects is 
due to the choice of the coupling which is of the same order as the on-site 
term (see the analysis in Ref. [20]). 

The two outermost rotobreathers systematically act as barriers between a cen- 
tral "hot" region of the chain and its boundaries, thus blocking the energy flow 
toward the environment. At some stage, such dynamical barriers are sponta- 
neously destabilized, thus allowing a portion of the trapped energy to rapidly 
flow away. This process goes on progressively reducing the size of the central 
hot patch until a single rotobreather survives in the bulk. We checked that its 
motion is periodic by computing the spatiotemporal spectrum S(k, uj) of the 
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Fig. 1. Relaxation in a chain with N = 200, 7 = 0.1, K = 1.4911, E(0)/N = 5. 
Upper panels G = 0, lower panels G = 1. (a), (c) Space-time density plots of the 
symmetrized site energies. The regions with higher energies are labeled with darker 
shades of grey. The chain sites i are reported on the horizontal axis (2 < i < N— 1), 
the vertical axis describes the time evolution in units of sampling times (250 natural 
time units), (b), (d) Decay of the normalized lattice energies. The time is expressed 
in units of the sampling time of the lattice site energies. 

(j)iS for several wave-numbers k. Indeed S(k, u>) displays sharp peaks at the 
fundamental frequencies plus small-amplitude harmonics. 

The right panels of Fig. 1 illustrate how the energy flows toward the reservoirs. 
In accordance with the above described phenomenology, the energy decay 
proceeds in a characteristic step-wise fashion, through a series of sudden jumps 
separated by long, approximately constant plateaus. The jumps correspond to 
the disappearance of one of the rotobreathers at the boundaries of the hot 
region. Between two consecutive jumps, the energy remains approximately 
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Fig. 2. Decay of the normalized total energy for four different initial conditions with 
E(0)/N = 5. (a) NN model, (b) Full sine-lattice model. The lattice size is N = 200. 
Other parameters are 7 = 0.1, K = 1.4911. 

constant with very small fluctuations. 

Our simulations show that this kind of scenario persists also in the case G = 0, 
thus confirming previous qualitative results [23]. This finding is also illustrated 
in Fig. 2, where we compare the energy decay of four different initial conditions 
for the two kinds of potentials. In order to emphasize the slow character of 
the process, we rescale the time by Tq = iV/27, which has been shown to be 
the characteristic time scale over which edge dissipation occurs [13]. 

It is interesting to observe that different initial conditions may give rise to 
substantially different sequences of jumps and plateaus, that may be regarded 
as different realizations of a stochastic process. For comparison, we remark 
that such sensitivity to the choice of the initial condition has not been ob- 
served in the relaxation of other non-linear systems, where different stories 
approximately resemble each other and display a much smoother decay [13,14]. 

In order to give a more quantitative description of the relaxation process, it 
is very interesting to study the average behavior of the energy decay. Fig- 
ure 3 shows that the latter is characterized by three different scaling regimes. 
These are best visualized by plotting the quantity — ln(E(t)/E(0)) in log- 
log scale, where a law of the form exp(— t a ) shows up as a straight line with 
slope a. We find that both the systems with G = and G = 1 exhibit the 
same properties. At short times (t <C r ), the energy decays exponentially as 
E(t) = E(0) exp(— t/r ). This is what one observes for a harmonic lattice of 
the same size, thus signaling that the non-linear effects are negligible on such 
a time scale. Actually, the same holds for the Fermi-Pasta-Ulam chain [13]. 
At intermediate times (t ~ 10~ 2 t ), the curves undergo a change of slope 
and slow down to stretched exponential laws. The knee associated with such 
crossover may be due to power-law corrections, which are known to be present 



6 






(b) 



Fig. 3. Decay of the normalized total energy averaged over a set of 99 initial con- 
ditions. We plot —\n(E(t)/E(0)) in log-log scale with E(0)/N = 5 and for two 
different lengths of the chain. The thin lines correspond to exp[— t/ro] (leftmost 
curve), and to a stretched exponential fit exp[— (i/r) "] of the data between the 
two arrows, (a) NN model, a = 0.59. (b) Full sine-lattice model a = 0.80. Other 
parameters are 7 = 0.1, K = 1.4911. 

at this stage in other systems [13]. Eventually (t > IOto), the residual quasi- 
stationary state is reached. In this last stage, a finite fraction of the initial 
energy is stored in a single rotobreather and remains constant over the rest of 
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Table 1 

Dependence of the exponent a on the initial energy density for N = 

E(0)/N G = G = 1 



100 and 7 = 0.1. 



3 



0.88 



0.93 



5 



0.59 



0.80 



7 



0.41 



0.63 



8 



0.40 



0.49 



the simulation. 

In Figure 3 we also consider the effects of the lattice size for the same value 
of the damping constant 7. In order to make the comparison, we rescale the 
time axis by the size-dependent constant tq = N/2 / y. We see that the data 
for N = 100, 200 collapse fairly well for t < r and follow approximately the 
same trend at later times. More marked differences are observed in the first 
crossover region. In particular, the range of stretched exponential decay seems 
to be larger for the longer chains. 

To complete this analysis, we checked the dependence of the exponent a on the 
initial energy. The data reported in Table 1 show that a decreases by increasing 
E(0)/N. Accordingly, this suggests that the higher the initial energy density, 
the longer is the "lifetime" of the rotobreathers that appear in the system. 
There are some results in the literature that support this inference for the 
sine-lattice systems (with and without local coupling). In fact, it has been 
found that the lifetime of an excitation of energy £ in a nonlinear system 
should be proportional to exp y/E [24]. 

Finally, we wish to mention that a similar stepwise decay process yielding 
an average stretched-exponential law has been found also in two-dimensional 
easy-plane ferromagnets [25]. In such case, annihilation of vortex-antivortex 
pairs is the basic physical mechanism. 



4 Statistical analysis 

In this section we attempt to understand better the origin of the slow kinet- 
ics. In particular, we shall focus on the statistics of plateaus and energy jumps 
during energy relaxation. As we already noticed above, the total energy of the 
lattice decreases in a stepwise fashion. Between subsequent jumps no apprecia- 
ble energy exchange with the external world occurs and a plateau is observed 
(see again Fig. 2), whose duration may vary over different time scales. For 
a single realization of a given initial condition, we can thus approximate the 
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curve E(t) as a piece-wise function through the sequences AE n (energy jump 
at the n-th step) and At n (duration of the n-th plateau): 



E(t) 



N p / n -l 



) 



E(0) 



l-J2AE n elt - ^ At 



n=l \ m=0 



rn 



(3) 



where Q(x) is the Heaviside step function. In the spirit of a statistical descrip- 
tion we thus regard each sequence (AE n , At n ) for n = 1, ... N p as a realization 
of a random process. 

Since all the features of this process are unknown, we devised a numerical 
procedure to extract from the simulations both the At n and AE n for each 
realization. The first task is to locate the beginning and end of each plateau to 
a reasonable level of accuracy. We chose to sample the whole time-series (with 
an adjustable sampling interval t s ) and look for the times where the energy 
differences J(t) = E(t + t s ) — E(t) exceed some pre-assigned threshold value 
\J(t)\ > q. Physically, J(t) is the energy flux toward the external reservoirs 
integrated over a time t s . 

In order to fix the value of the threshold, we computed the average number 
of plateaus (N p ) as a function of q. We found that there is always a range of 
q values where (N p ) displays a broad maximum. Therefore, we systematically 
chose the value of q in such a way as to minimize the variation of the number of 
steps. According to this prescription, our procedure identified approximately 
the same number of relaxation stages at each realization. The average value of 
the number of plateaus was found to depend primarily on the type of potential 
energy of the lattice, while displaying little sensitivity on the initial energy 
density. For example, with N = 160 the average value of the relaxation steps 
before reaching the quasi-stationary state in the NN lattices is (N p ) fa 11 ± 1, 
while almost twice as much relaxation stages are needed in the presence of on- 
site term, (N p ) fa 26±1. We infer that the latter system could be characterized 
by a much rougher energy landscape. 

Since we do not know whether the process at hand is a stationary one, we 
first checked how the averages depend on time, i.e. on the index n. In Fig. 4 
we see that the averages (At n ) are almost independent of the order of the 
relaxation n for both G = and G ^ 0. At least to a first approximation, it is 
thus legitimate to assume that the process is a stationary one. It is important 
to stress that the above observation rules out the conjecture that the slow 
decay kinetics of the system may be the result of a progressive increase of the 
residence times in metastable states visited during relaxation. 

The next step is to investigate the probability distributions of the lifetimes 
V(At) and energy jumps V(AE). In Fig. 5 (a) we show the histograms of the 
energy jumps for the two cases G = and G^O. We found that the data can 
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Fig. 4. Plot of the average duration of the plateaus (At n ) vs order of the relaxation 
stage n for G = and G = 1. Parameters are: E(0)/N = 8, N = 160, 7 = 0.1. 

be fitted by the empirical normalized function 



Hence, the distribution of V(AE) is exponentially small for small values of AE 
and decays following a power law of the kind l/AE a+1 for large values of the 
energy jumps. In the case shown in Fig. 5, we found a ~ 1.6, AE 0.0081 
(G = 0), and a « 1.8, AE « 0.0074 (G = 1). 

We turn now to the distribution of plateaus V(At). In practice, it is a difficult 
task to construct the histograms from the measured time series, since the latter 
span a wide range (3 decades) yielding exceedingly large fluctuations. A way 
around this problem is to construct the cumulative (integrated) distributions 
C(At), that can be obtained from the numerical series by inverting the axes in 
a simple rank-size plot of the data. Remarkably, we found that the statistics of 
the time intervals follows an inherently slow stretched exponential distribution 
law. This finding is illustrated in Fig. 5 (b) for both the G = and G ^ 
cases, where we plot the quantity — log[C(At)]. More precisely, we see that the 
first part of the distribution is a pure exponential (straight lines with slope one 
in Fig. 5 (b)), while the tails follow a stretched exponential law. This means 
that the probability distribution V(At) (the derivative of C(At)) is 



up to power-law corrections. The stretched exponential portion of the curves 
sets in at times longer than At* 10 3 , which is systematically below the 
average, approximately constant duration of each relaxation stage. This ob- 
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Fig. 5. Statistics of the relaxation process, (a) Normalized histograms of energy 
jumps AE for G = (triangles) and G = 1 (circles) and two-parameter fit with 
the expression (4). (b) Cumulative distribution of the plateau durations for G = 
and G = 1 (symbols). The dashed lines are pure exponential fits of the portion 
At < At*. The solid lines are stretched exponential fits in the region At > At* 
(slightly shifted for the sake of clarity). The best -fit values of the exponent (5 are: 
f3 = 0.43 (G = 0) and (3 = 0.52 (G = 1). The series {At}, {AE} have been collected 
on a total of 99 initial conditions. Parameters are: E(0)/N = 8, N = 160, threshold 
q = 0.001. 
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servation suggests that the statistics of the plateau duration is dominated by 
the stretched exponential tail. Such intrinsically fat-tailed distribution of a 
series of relaxation events is thus directly related to the slow energy decay. 

In general, the origin of a stretched exponential relaxation is traced back to 
a superposition of many different time scales, distributed according to given 
weights [12,26]. In our context, however, it is not clear how to identify such 
a hierarchical set of degrees of freedom, that would be responsible for the 
observed slow kinetics. It is thus legitimate to explore the possibility of alter- 
native explanations. An interesting perspective has been presented by Bunde 
et al. in a recent paper [27], in connection with the statistics of extreme 
events. According to the authors, a stretched exponential law may rule the 
distribution of intervals between consecutive rare events provided the origi- 
nal time series displays power-law correlation. More precisely, they consider 
a stochastic time series x(t), whose autocorrelation decays by construction as 
t~ 7 with < 7 < 1. In this context, rare events are defined by the condition 
x(t) > q, where q is a pre-assigned threshold. Their numerical analysis shows 
that the distribution of the time intervals At separating extreme events decays 
as exp[— At 7 '], with 7' ss 7. 

It appears natural to extend this kind of description to the process of energy 
release in our system: each energy step AE, marking the collapse of a roto- 
breather, can be regarded clS cl TdbTG event " (sharp peak) in the series of energy 
flux J {t) (see Fig. 2). In order to validate this analogy, we checked that long- 
range correlations do exist in the series J(t) (Fig. 6). Within the statistical 
fluctuations, the autocorrelation function of J(t) decreases approximately as 
an inverse power law, whose exponent is smaller than one (~ 0.7). This value 
is of the same order of the exponents reported in Fig. 5 (b). As a matter of 
fact, the statistical fluctuations prevent us from obtaining a better estimate. 
Furthermore, the distribution of plateau durations are only asymptotically 
stretched-exponential, thus introducing a further source of error. Therefore, 
to make our analogy more quantitative further numerical work should be car- 
ried out. 

We note that a similar analysis of the correlations in the series of the energy 
flux in the Fermi-Pasta-Ulam chain does not yield a power-law decay. In fact, 
the autocorrelation function of J(t) vanishes exponentially. This is consistent 
with the absence of stretched exponential relaxation in this model [13], thus 
suggesting that the microscopic mechanisms operating in the chain of coupled 
rotators are substantially different. 

In order to provide further support to our interpretation, we have set up a 
Monte Carlo procedure which generates stochastic series of N p energy jumps 
and plateaus as in Eq. (3) according to given distributions V(AE) and V(At). 
We have then computed the average over such a set of realizations according 



12 




10 



10 10 
time lag 



10 



Fig. 6. Autocorrelation of the time series J(t) averaged over 99 initial conditions 
for N = 160, 7 = 0.1, E(0)/N = 8. The dashed line is an inverse power law with 
exponent 0.7. 

to Eq. (3), where N p can be fixed close to its empirical average value. We found 
that, by using a distribution of the form (4) and a stretched exponential law 
for the V(At), we could recover a stretched exponential decay of the average 
energy by just tuning at most two of the free fitting parameters AEq, a, T2, j3. 
For example, by fixing N p = 15, AE = 0.0081 and a = 1.7 in the case 
G — 0, E(0)/N = 8, the best-fit value of the exponent a from the Monte 
Carlo simulation [a = 0.42) is found to be in very good agreement with the 
one directly fitted on the simulation data (see Table 1). 



5 The effect of isotopic disorder 



In this section we investigate the role of disorder on the relaxation dynamics of 
the sine-lattice system. More specifically, we consider model (1) with randomly 
distributed momenta of inertia, that is either 7j = R (with probability p) or 
— 1 (with probability l—p). This choice corresponds to a 4-methyl-pyridine 
quasi-lD chain with random insertions of deuterated methyl groups, that has 
been studied experimentally by Fillaux and co-workers [16,28]. Hence, we 
fix R = 2, that corresponds to the ratio of momenta of inertia reported in 
Ref. [16]. 

The motivation of this analysis is twofold. First, the joint effect of non-linearity 
and disorder on the non-equilibrium properties of discrete lattices are far 
from being satisfactorily understood [29]. Second, the disordered sine-lattice 
model has been used to interpret data from inelastic neutron scattering ex- 
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Fig. 7. Decay of the normalized total energy in a chain with isotopic disorder. We 
plot -ln(E(t)/E(0)} in log-log scale with E(0)/N = 4.6 and for three different 
values of the concentration p. The data have been averaged over a set of 30 initial 
conditions. Parameters are N = 100, G = 1, 7 = 0.1. 

periments [16]. More specifically, it has been shown that the effects of isotopic 
disorder on the dynamics rapidly disappear upon increasing the temperature. 
This observation has been explained quantum-mechanically as a transition 
to the classical regime, thus suggesting that the effects of disorder should be 
negligible in the latter case. By the same token, one may expect that the re- 
laxation dynamics in the classical system should display weak sensitivity on 
the concentration of impurities. 

Our simulations confirm the above hypothesis. In Fig. 7 we show that the en- 
ergy decay is hardly affected by the value of p. Remarkably, the curves almost 
overlap and display both a wide region of stretched exponential behavior and 
the saturation to a residual state. Accordingly, the estimated exponents are 
fairly insensitive to the concentration of impurities. Moreover, we checked that 
the energy which is left in the residual state is approximately independent of 
p, with deviations of at most 10 % from the average. 



6 Conclusions 

The study of the energy relaxation in models of classical rotators has allowed 
us to identify a possible mechanism yielding the stretched exponential behav- 
ior, observed also in other one-dimensional models [11,12,29]. The static na- 
ture of localized breather states, which are spontaneously formed upon cooling 
the lattice from its boundaries, is one of the basic ingredients. Indeed, static 
localized solutions can form "dynamical barriers" that are almost decoupled 
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from their neighbours and thus segregate an internal hot lattice region from the 
cooled boundaries. The energy initially fed into the system remains trapped 
there, until the rotobreathers are abruptly destroyed. It may be argued that 
this sudden death must originate from a sufficiently large resonant fluctuation 
emerging from the hot core. In our approach, we have employed a statistical 
description, which does not make any explicit reference to specific dynamical 
mechanisms. 

Accordingly, we showed that the destruction of boundary rotobreathers may 
be regarded as a "rare event". The resulting statistics depends on two main 
factors: the energy of the rotobreathers and the probability that in a given time 
interval the hot region may produce an excitation capable of annihilating them. 
In our numerical analysis we can just observe the effects of this complicated 
interaction, which results in a stretched-exponential law for the energy decay. 

It may be argued that the mechanism of energy trapping must be rather 
peculiar to the ID nature of the model. Indeed, we expect that in higher- 
dimensional lattices this phenomenon may not occur at all, and that the char- 
acteristic stepwise decay of energy could be replaced by a smoother decrease. 
Some preliminary simulations of a 2D arrays of rotators confirm this idea. 
In particular, we observe that the residual state resembles the one of the 2D 
Fermi-Pasta-Ulam model - a collection of static breathers randomly arranged 
on the lattice [14]. 

Altogether, we have seen that the possibility of easily exciting long-lived roto- 
breathers is responsible for the slow kinetics. It should be recalled that these 
very same objects have been invoked to explain the normal conduction in the 
NN model [19]. In this sense, our results complement those observations and 
confirm the essential role of nonlinear excitations in the out-of-equilibrium 
properties of many-body systems. 
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